%% Carbon Taxes Around the World: Cooperation, Strategic Interactions, and Spillovers
% IMF Economic Review
% Alessandro Moro and Valerio Nispi Landi
% Replication files
% This file computes the reaction functions and the Nash equilibrium

clc; clear all; close all;
RUN=1;


if RUN==1 % The Dynare loop starts to find welfare for each policy combination

tic
console
console_final_opt
if flag<1
disp(['Computation fails in console']);
return
end
BASE=0;
SAVE=1;
dynare climate



%% Find the Nash equilbrium

[~,reaction_ind1] = max(WW);        % maximum of the columns of WW
[~,reaction_ind2] = max((WWz'));    % maximum of the rows of WWz
reaction1=DG(reaction_ind1);
reaction2=DG(reaction_ind2);
reaction=[reaction1;reaction2;DG]';


ne1=reaction(1,1);
ne2=reaction(1,2);
limit=1000;
j=1;

while limit>0
ind1=find(DG==ne2(j));
ne1(j+1)=reaction(ind1,1);
ind2=find(DG==ne1(j+1));
ne2(j+1)=reaction(ind2,2);
limit1=abs(ne1(j+1)-ne1(j));
limit2=abs(ne2(j+1)-ne2(j));
limit=norm([limit1;limit2]);
j=j+1;
end

ne=[ne1(end),ne2(end)];
% Impose symmetry
if SYM==1 && ind1~=ind2
ind2=ind1;
ne=[ne1(end),ne1(end)];
end

WW_ne_max=WW(ind1,ind2);
WWz_ne_max=WWz(ind1,ind2);

%% Find the Coordination equilbrium

WWw_bis=WWw;

% Add the constraint that countries are not worse off compared to the Nash equilibrium
AE=WW>WW_ne_max;
EME=WWz>WWz_ne_max;

for rr=1:length(DG)
for cc=1:length(DGZ)
if AE(rr,cc)==0 || EME(rr,cc)==0
WWw_bis(rr,cc)=-1e9;
end
end
end

[~,ind_w] = max(WWw_bis(:));
[dd,ee] = ind2sub(size(WWw_bis),ind_w); % find the position of the optimal coefficients under the coordination equilibrium

% Impose symmetry
if SYM==1 && dd~=ee
ee=dd;
end
coord(1)=DG(dd);   % Optimal tau
coord(2)=DGZ(ee);  % Optimal tauz

WW_coord_max=WW(dd,ee);
WWz_coord_max=WWz(dd,ee);

%% Find the social-planner equilbrium


[~,ind_w] = max(WWw(:));
[dd,ee] = ind2sub(size(WWw),ind_w); % find the position of the optimal coefficients under the coordination equilibrium

% Impose symmetry
if SYM==1 && dd~=ee
ee=dd;
end

sp(1)=DG(dd);   % Optimal tau
sp(2)=DGZ(ee);  % Optimal tauz
WW_sp_max=WW(dd,ee);
WWz_sp_max=WWz(dd,ee);

cequiv=exp((1-beta)*(WW_coord_max-WW_ne_max))-1;
cequivz=exp((1-betaz)*(WWz_coord_max-WWz_ne_max))-1;

cequiv_sp=exp((1-beta)*(WW_sp_max-WW_ne_max))-1;
cequivz_sp=exp((1-betaz)*(WWz_sp_max-WWz_ne_max))-1;
%% Saving ne and coord

if psiX==0 && SAVE==1
WW_psi0=WW;
WWz_psi0=WWz;
WWw_psi0=WWw;
cequiv_psi0=cequiv;
cequivz_psi0=cequivz;
cequiv_sp_psi0=cequiv_sp;
cequivz_sp_psi0=cequivz_sp;
reaction1_psi0=reaction1;
reaction2_psi0=reaction2;
ne_psi0=ne;
coord_psi0=coord;
sp_psi0=sp;

save sim_psi0 WW_psi0 WWz_psi0 WWw_psi0 reaction1_psi0 reaction2_psi0 ne_psi0 coord_psi0...
      cequiv_psi0 cequivz_psi0 cequiv_sp_psi0 cequivz_sp_psi0 sp_psi0
elseif psiX==0.5 && SAVE==1
WW_psi05=WW;
WWz_psi05=WWz;
WWw_psi05=WWw;
cequiv_psi05=cequiv;
cequivz_psi05=cequivz;
reaction1_psi05=reaction1;
reaction2_psi05=reaction2;
ne_psi05=ne;
sp_psi05=sp;
cequiv_sp_psi05=cequiv_sp;
cequivz_sp_psi05=cequivz_sp;
coord_psi05=coord;

save sim_psi05 WW_psi05 WWz_psi05 WWw_psi05 reaction1_psi05 reaction2_psi05 ne_psi05 coord_psi05...
     cequiv_psi05 cequivz_psi05 cequiv_sp_psi05 cequivz_sp_psi05 sp_psi05
elseif psiX==1 && SAVE==1
WW_psi1=WW;
WWz_psi1=WWz;
WWw_psi1=WWw;
cequiv_psi1=cequiv;
cequivz_psi1=cequivz;
reaction1_psi1=reaction1;
reaction2_psi1=reaction2;
ne_psi1=ne;
sp_psi1=sp;
cequiv_sp_psi1=cequiv_sp;
cequivz_sp_psi1=cequivz_sp;
coord_psi1=coord;

save sim_psi1 WW_psi1 WWz_psi1 WWw_psi1 reaction1_psi1 reaction2_psi1 ne_psi1 coord_psi1...
    cequiv_psi1 cequivz_psi1  cequiv_sp_psi1 cequivz_sp_psi1 sp_psi1
end

% Load the sim_psi0/psi1/psi05 files

elseif RUN==0 
console
SAVE=0;
LIN=0;
DAM=0;
load_sim
end


disp(['-----------------------------------------------------------------------']);
disp(['Nash eq: AE: Ϝ=',num2str(ne(1)*100),   '%,  Ϝ*=',num2str(ne(2)*100),'%']);
disp(['SP eq: Ϝ=' ,num2str(sp(1)*100),'%, Equiv ',num2str(cequiv_sp*100) '%; Ϝ*=',num2str(sp(2)*100),'%, Equiv ',num2str(cequivz_sp*100) '%']);
disp(['Coor eq: Ϝ=' ,num2str(coord(1)*100),'%, Equiv ',num2str(cequiv*100),'%; Ϝ*=',num2str(coord(2)*100),'%, Equiv ',num2str(cequivz*100) '%']);
disp(['-----------------------------------------------------------------------'])

